EXACT PHYLODYNAMIC LIKELIHOOD VIA STRUCTURED MARKOV GENEALOGY PROCESSES

We consider genealogies arising from a Markov population process in which individuals are categorized into a discrete collection of compartments, with the requirement that individuals within the same compartment are statistically exchangeable. When equipped with a sampling process, each such population process induces a time-evolving tree-valued process defined as the genealogy of all sampled individuals. We provide a construction of this genealogy process and derive exact expressions for the likelihood of an observed genealogy in terms of filter equations. These filter equations can be numerically solved using standard Monte Carlo integration methods. Thus, we obtain statistically efficient likelihood-based inference for essentially arbitrary compartment models based on an observed genealogy of individuals sampled from the population.


Introduction.
When the genome of an infectious agent accumulates mutations on timescales similar to those of transmission and infection progression, the resulting pattern of differences among genomes contains information on the history of the pathogen's passage through individual hosts and the host population.As Grenfell et al. (2004) observed, one can extract this information to gain insight into the structure and dynamics of the host-pathogen system.In particular, one can formalize mathematical models of transmission, estimate their parameters, and compare their ability to explain data, following standard statistical paradigms.This is known as phylodynamic inference.Alizon (2024) gives a good review of the history of the subject.
The most common approach to phylodynamic inference rests upon a mathematical linkage between the tree-like genealogy or phylogeny that expresses the relationships of shared ancestry among sampled genomes and a model of the dynamics of the transmission system.Various linkages are possible, but because it is maximally efficient (i.e., loses the least information), it is desirable to be able to compute the likelihood function for models of interest.This is simply the probability density of a given genealogy conditional on a given model, viewed as a function of the parameters of that model.In particular, if S is a set of genome sequences, Φ a genealogical tree relating these sequences, E a model of sequence evolution, and D a dynamic transmission model, then the likelihood is LpD, Eq " f pS|D, Eq " ż f pS|Φ, Eq f pΦ|Dq dΦ, where the integral is taken over all possible genealogies and we somewhat loosely use the symbol f for the various distinct probability densities, the nature of each of which is clear from its arguments.In this expression, f pS|Φ, Eq is typically the Felsenstein (2004) phylogenetic likelihood.The function f pΦ|Dq, which links the phylogeny to the dynamic model, may be termed the phylodynamic likelihood.In the Bayesian context, this same function is sometimes referred to as a tree prior (Möller et al., 2018;Volz & Siveroni, 2018).The computation of this function has remained out of reach, except in several special cases.This paper presents theory that enables its computation for a very broad range of dynamic models.
Existing approaches to the phylodynamic likelihood have been based on one of two mathematical idealizations.
The first is the Kingman (1982a) coalescent, by which likelihood of a given genealogy is computed using a reverse-time argument.This computation provides the exact likelihood for a genealogy resulting from a particular, constant population-size, dynamic model (the Moran model, e.g., Moran, 1958;Kingman, 1982b;Möhle, 2000).Extensions of this approach develop approximate likelihoods for the case when the population size varies as a function of time (Griffiths & Tavaré, 1994;Drummond et al., 2005) or according to an SIR process (Volz et al., Date: May 28, 2024. 2009; Rasmussen et al., 2011), as long as the population size is large and the sample-fraction remains negligible.
The second idealization is the linear birth-death process, for which exact expressions for the likelihood are available (Stadler, 2010).Linearity in this context amounts to the assumption that distinct lineages do not interact: it is the resulting self-similarity of genealogies that renders the likelihood analytically tractable.Extensions of this approach develop approximations via linearization of nonlinear processes or restriction to scenarios in which population growth is nearly linear (e.g., MacPherson et al., 2021).Although the tractability of these approaches makes them attractive, concern naturally arises as to validity of the approximations in specific cases, the biases introduced by them, and the amount of information in data left unutilized by these approximate methods.For this reason, there is interest in improved phylodynamic inference techniques.
What would an ideal phylodynamic inference method look like?First, it would afford exact computation of the phylodynamic likelihood, so that comparisons among parameterizations and models could be made on a sound basis.Second, because nonlinearity, nonstationarity, noise, and measurement error are prominent and ubiquitous in epidemiology, it would accommodate nonlinear, time-inhomogeneous, stochastic transmission models.Third, because many of the most scientifically important uncertainties concern heterogeneities in transmission rates and the susceptibility, behavior, age, and location of hosts, it would accommodate host populations structured by these factors.While some structuring factors (e.g., age, spatial location) are most naturally expressed in terms of continuous variables, discretely structured models have repeatedly proved their value in epidemiology.
In particular, compartmental models are extremely flexible and have often been used as approximations when continuous structure leads to uncomfortably high model dimension.Finally, because there is typically uncertainty not only in the parameters, but also in the structure, of a host-pathogen system, an improved phylodynamic inference methodology would place minimal restrictions on the form of the models that it can accommodate.This paper demonstrates how these desiderata can be achieved-at least for models with discrete structure-including arbitrary nonlinear compartmental models with a countable number of compartments.
To connect a model at the level of a population with genealogies based on samples taken from individual hosts, it is necessary to make assumptions about the individuals in the population.The simplest such assumption is that the individuals that are identical with respect to the population dynamics are indeed statistically identical.That is, that they are exchangeable.In a compartmental model, this is tantamount to the assumption that the residence times of the individuals within each compartment are identically distributed, though not independent.Although exchangeability is indeed an additional assumption, it is so natural that it is frequently unrecognized as such, and one often reads statements to the effect that exchangeability of individuals is a consequence of the Markovian assumption.Nonetheless, since it adds minimal additional structure, it is the natural assumption, and the one we will make in this paper.
the left-limit of such a process.Accordingly, if Φ t is a càdlàg random process, we define Note that r Φ t is thus left-continuous with right limits.
If Φ t , t P R `is a pure jump process, knowledge of its sample path is equivalent to knowledge of the number, K t , of jumps it has taken as of time t, the jump times Tk , and the embedded chain Φk -Φ Tk , k " 0, . . ., K t .
2.2.Population process.We are motivated by the desire for exact phylodynamic inference methods for as wide a class of epidemiological models as possible.In particular, we would like to be able to formulate and parameterize an arbitrary compartmental model and to quantify its ability to explain data using likelihood.Fig. 1 depicts a few such models in order to give a sense of the kinds of complexities that can arise.Of course, with the ability to entertain models with countably many compartments, much greater complexity is possible.In particular, one can model not only complex infection progression, but also strain structure, behavioral structure, age structure, and spatial structure using compartmental models.As is well known, one can discretize continuous structurevariables and employ the linear chain trick to accommodate non-exponential residence times.While the utility of these approximations will vary, a very wide range of model assumptions lie within the scope of the theory presented here.
We will assume that our population process is a time-inhomogeneous Markov jump process, X t , t P R `, taking values in some space X.In earlier work (King et al., 2022), we limited ourselves to the case X " Z d , but here we assume only that X is a complete metric measure space with a countable dense subset.The population process is completely specified by its initial-state density, p 0 , and its transition rates α.In particular, we suppose that (1) Prob rX 0 P Es " ż E p 0 pxq dx for all measurable sets E Ď X.For any t P R `, x, x 1 P X, we think of the quantity αpt, x, x 1 q as the instantaneous hazard of a jump from x to x 1 .More precisely, the transition rates have the following properties: αpt, x, x 1 q ě 0, ż X αpt, x, x 1 q dx 1 ă 8, for all t P R `and x, x 1 P X and that, as a function of time, α is continuous almost everywhere.Henceforth, we understand that integrals are taken over all of X unless otherwise specified.Let K t be the number of jumps that X has taken by time t.We assume that K t is a simple counting process so that We further assume that αpt, x, x 1 q is càdlàg as a function of time for all x, x 1 P X and that the number of jumps that occur in a finite time-interval is finite, i.e., Prob rK t ă 8s " 1 for all t.
2.3.Kolmogorov forward equation.The above may be compactly summarized by stating that if vpt, xq satisfies the Kolmogorov forward equation (KFE), (2) Bv Bt pt, xq " ż vpt, x 1 q αpt, x 1 , xq dx 1 ´ż vpt, xq αpt, x, x 1 q dx 1 , and if, moreover, vp0, xq " p 0 pxq, then ş E vpt, xq dx " Prob rX t P Es for every measurable E Ď X. Eq. 2 is sometimes called the master equation for X t .2.4.Inclusion of jumps at deterministic times.For modeling purposes, it is sometimes desirable to insist that certain events occur at known times.For example, if samples are collected at specific times in such a way that the timing itself conveys no information about the process, one might wish to condition on the sampling time.We can expand the class of population models to allow for this as follows.Suppose that S " ts 1 , s 2 , . . .,u Ă Z ìs a sequence of event times.Let us postulate that, at each of these times, an event occurs at which X t jumps according to a given probability kernel π.In particular, for any state x P X and measurable E Ă X, πps i , x, Eq is the probability that the jump at time s i is to E, conditional on the state just before the jump being x.With this notation, the KFE for the process becomes Note that the Eq. 3 is identical to Eq. 2; we call this the regular part of the KFE.We refer to Eq. 4 as the singular part of the KFE.
x " pS, E, I, Rq x 1 " pS ´1, E `1, I, Rq x 1 " pS, E ´1, I `1, Rq x 1 " pS, E, I ´1, R `1q x 1 " pS `1, E, I, R ´1q x 1 " pS, E, I, Rq As a matter of notation, one can represent Eqs. 3 and 4 as a single equation in the form of Eq. 2. In particular, if in Eq. 2 we make the substitution αpt, x, x 1 q Þ Ñ αpt, x, x 1 q `ÿ sPS δpt, sq dπ dx 1 pt, x, x 1 q, we obtain an equation which we can view as shorthand for Eqs. 3 and 4. Here, δpt, sq is a Dirac delta function and dπ{dx 1 denotes the density (i.e., Radon-Nikodym derivative) of π with respect to the measure on X.
2.5.Jump marks.It will be useful to divide the jumps of the population process X t into distinct categories, which differ with respect to the changes they induce in a genealogy.For this purpose, we let U be a countable set of jump marks such that αpt, x, x 1 q " ÿ uPU α u pt, x, x 1 q.
Fig. 2 shows an example for which U has five elements.In the following, sums over u are to be taken over the whole of U unless otherwise indicated.
Let us define the jump mark process, U t , to be the mark of the latest jump as of time t.As usual, we take the sample paths of U t to be càdlàg.Observe that, though X t and pX t , U t q are Markov processes, U t is not.
2.6.Demes and deme occupancy.Our first goal in this paper is to show how a given population process induces a unique stochastic process on the space of genealogies.At each time, this genealogy will represent the relationships of shared ancestry among a population of lineages extant at that time.To accommodate the structure of the population, this population of lineages will itself be subdivided into discrete categories.In particular, we suppose that there are a countable set of subpopulations, within each of which individual lineages are exchangeable.We call these subpopulations demes, and use the symbol D to denote an index set for them.Fig. 1 illustrates this concept in the context of several compartmental models.
We define the deme occupancy function n : D ˆX Ñ Z `so that for i P D, x P X, n i pxq is the number of lineages in deme i when the population is in state x.

Examples.
The class of population models to which the theory presented here applies is very broad indeed.
In particular, it encompasses the entire class of compartmental models with time-dependent flow rates.Here, to give a sense of this breadth, we briefly describe a few models of interest.Appendix B works out the theory for each of these examples.
SIRS model.King et al. (2022) worked out formulas for the exact likelihood of a genealogy induced by an SIRS model.The theory developed in this paper applies, but since there is only one deme in this model, this is a simple case.
SEIRS model.A simple, yet interesting, model with more than one deme is the SEIRS model (Fig. 1A).The state space is Z 4 `, with the state x " pS, E, I, Rq defined by the numbers of hosts in each of the four compartments.It has two demes: D " tE, Iu.The deme occupancy function in this case is npxq " pE, Iq.Note that the terms associated with sampling cancel each other in the KFE, since, in this model, sampling has no effect on the state.
Two-strain competition model.A simple model for the competition of two strains for susceptible hosts is depicted in Fig. 1B.In this model, the state vector consists of seven numbers: Superspreading model.Linear birth-death model.The linear birth-death process, a mainstay of existing phylodynamic methods, is a special case of the theory presented here.For this process, we have X " Z `and there is a single deme.X t represents the size of a population and npX t q " X t .
Moran model and the Kingman coalescent.The Kingman (1982a) coalescent is another workhorse in existing phylodynamic approaches.It is the ancestral process for the Moran model, in which a fixed population of n lineages experiences events at times distributed according to a rate-µ Poisson process.At each such event, an individual lineage selected uniformly at random dies and is replaced by the offspring of a second randomly selected lineage.
2.8.History.Consider the Markov process pX t , U t q.We define its history process, H t , to be the restriction of the random function s Þ Ñ pX s , U s q to the interval r0, ts.Note that H t is itself trivially a Markov process, since it contains its own history.
Alternatively, one can think of H t as consisting of the sequence ´´T k , Xk , Ûk ¯¯Kt k"0 . In particular, conditional on H t , both X t and U t are deterministic, as are K t , the embedded chains, Xk , Ûk , and the point process of event times Tk .The probability measure on the space of histories can be expressed in terms of these: (5) Prob rdHts " p0p X0q d X0 where again, by convention, T0 " 0 and TK t `1 " t.
If H is such a history, we define tpHq to be the right endpoint of its domain and use the notation evpHq -! T1 , . . ., TKt ) Ă r0, tpHqs to denote the set of its jump times.
2.9.Genealogies.A genealogy, G, encapsulates the relationships of shared ancestry among a set of lineages that are extant at some time tpGq P R `and perhaps a set of samples collected at earlier times (Fig. 3).A genealogy has a tree-or forest-like structure, with four distinct kinds of nodes: (i) tip nodes, which represent labeled extant lineages; (ii) internal nodes, which represent events at which lineages diverged and/or moved from one deme to another; (iii) sample nodes, which represent labeled samples; and (iv) root nodes, at the base of each tree.Each node a is associated with a specific time, tpaq.In particular, if a is a tip node in G, then tpaq " tpGq; if a is a sample node, then tpaq ď tpGq is the time at which the sample was taken.Moreover, if node a is ancestral to node a 1 , then tpaq ď tpa 1 q and tpa 1 q ´tpaq is the distance between a and a 1 along the genealogy.Without loss of generality we assume that tpaq " 0 for all root nodes a.We let evpGq denote the set of all internal and sample node-times of the genealogy G; we refer to these as genealogical event times.
Importantly, a genealogy informs us not only about the shared ancestry of any pair of lineages, but also about where in the set of demes any given lineage was at all times.Accordingly, we can visualize a genealogy as a tree, the nodes and edges of which are painted with a distinct color for each deme (Fig. 3).Note that a genealogy will in general have branch-point nodes, i.e., internal nodes with more than one descendant, but may also have internal 0.0 2.5 5.0 7.5 10.0 FIGURE 3. A genealogy, G, specifies the relationships of shared ancestry (via its treestructure) and deme occupancy histories (via the coloring of its branches) of a set of lineages extant at some time tpGq, as well as some samples gathered at earlier times.Here, tpGq " 10 and there are two demes, D " tblue, yellowu.Tip nodes, denoting extant lineages, are shown as black dots; sample nodes are shown as blue dots; internal nodes are indicated in green.Note that internal nodes occur not only at branch-points, but also inline (i.e., along branches).Wherever a lineage moves from one deme (color) to another, an internal node occurs; the converse does not necessarily hold.
nodes with only one descendant.We refer to such nodes as inline nodes.These occur whenever the color changes along a branch, but can also occur without a color-change.
Formally, we define a genealogy, G, to be a triple, pT, Z, Y q, where T " tpGq P R `is the genealogy time, Z specifies the genealogy's tree structure, and Y gives the coloring.In particular, let L be a countable set of labels and let partitpLq be the set of all collections of finite, mutually-disjoint subsets of L. That is, an element z P partitpLq is a partition of the finite set Ť z Ď L. Partition fineness defines a partial order on partitpLq.Specifically, for z, z 1 P partitpLq, we say z ≼ z 1 if and only if for every b 1 P z 1 there is b P z such that b Ě b 1 .The tree structure of G is defined by a càdlàg map Z : r0, T s Ñ partitpLq that is monotone in the sense that t 1 ď t 2 implies Z t1 ≼ Z t2 .An element b P Z t is a set of labels; it represents the branch of the tree that bears the corresponding lineages.We use the notation evpZq to denote the set of times at which Z is discontinuous.Note that evpZq includes the times of all tip, sample, and branch-point nodes, but excludes inline and root nodes.Therefore, evpZq Ď evpGq.
The third element of G specifies the coloring of branches and locations of tip, sample, and internal nodes (including inline nodes).Mathematically, if G " pT, Z, Y q, then Y is a càdlàg function that maps each point on the genealogy to a deme and a non-negative integer.In particular, if t P r0, T s and a is the label of any tip or sample node, Y t paq " pY d t paq, Y m t paqq P D ˆZ`, where Y d t paq is the deme in which the lineage of a is located at time t and Y m t paq is the number of internal or sample nodes encountered along the lineage of a in going from time 0 to time t.In particular, Y m t paq is a simple counting process, with Y m 0 paq " 0 for all a.Since a, a 1 P b P Z t implies Y t paq " Y t pa 1 q, one can equally well think of Y t as a map Z t Ñ D ˆZ`.Given a tree Z, we let YpZq denote the set of colorings Y that are compatible with Z.We moreover define Y t pZq -tY t | Y P YpZqu.Formally speaking, YpZq is a fiber bundle over Z, each Y t pZq being a fiber.
It will sometimes be convenient to make use of notation whereby a genealogy G " ptpGq, G Z , G Y q.
2.10.Binomial ratio.For n, r, ℓ, s P Z D `, define the binomial ratio Observe that ¨n ℓ r s ‹ ‹ ‚ P r0, 1s.Moreover, in consequence of the Chu-Vandermonde identity, we have 3. The induced genealogy process.
3.1.Event types.We now show how a given population process naturally induces a process in the space of genealogies.Specifically, at each jump in the population process, a corresponding change occurs in the genealogy, according to whether lineages branch, die, move between demes, or are sampled.For this purpose, there are five distinct pure types of events: (a) Birth-type events result in the branching of one or more new lineages, each from some existing lineage.
Examples of birth-type events include transmission events, speciations, and actual births.Importantly, we assume that all new lineages arising from a birth event share the same parent and that at most one birth event occurs at a time, almost surely.(b) Death-type events result in the extinction of one or more lineages.Examples include recovery from infection, death of a host, and species extinctions.We allow for the possibility that multiple lineages die simultaneously.(c) Migration-type events result in the movement of a lineage from one deme to another.Spatial movements, changes in host age or behavior, and progression of an infection can all be represented as migration-type events.We permit multiple lineages to move simultaneously.(d) Sample-type events result in the collection of a sample from a lineage.We allow for the possibility that multiple samples are collected simultaneously, though we require that, in this case, each extant lineage is sampled at most once.(e) Neutral-type events result in no change to any of the lineages.Fig. 2 depicts an example with jumps of all five pure types.It is not necessary that an event be of a pure type; compound events partake of more than one type.For example, a sample/death-type event, in which a lineage is simultaneously sampled and removed, has been employed (Leventhal et al., 2014), as have birth/death events in which one lineage reproduces at the same moment that another dies (e.g., the Moran (1958) process).The theory presented here places few restrictions on the complexity of the events that can occur by combining events of the various pure types.
3.2.Genealogy process.We now show how a given population process induces a stochastic process, G t , on the space of genealogies.In the case of unstructured population processes (i.e., those having a single deme), King et al. (2022) gave a related construction that is equivalent to the one presented here.
At each jump in the population process, a change is made to the genealogy, according to the mark, u, of the jump (Fig. 4).In particular: (a) If u is of birth-type (Fig. 4A), it results in the creation of one new internal node, call it b.A tip node, a, of the appropriate deme is chosen with uniform probability from among those present and b is inserted so that its ancestor is that of a, while a takes b as its ancestor.One new tip node, of the appropriate deme, is created for each of the children, all of which take b as their immediate ancestor.(b) If u is of death-type (Fig. 4B), one or more tip nodes of the appropriate demes are selected with uniform probability from among those present.These are deleted.Next, internal nodes without children are recursively removed.Sample nodes are never removed.(c) At a migration-type event (Fig. 4C), the appropriate number of migrating lineages are selected at random with uniform probability, from among those present in the appropriate demes.For each selected lineage, one new branch node is inserted between the selected tip node and its ancestor.The color of the descendant branch changes accordingly.
Event types differ by their effects on the genealogy.This can be seen by examining the local structure of the genealogy in the neighborhood of a jump.(A) A birth-type jump results in the branching of one or more child lineages from the parent.There can be only one parent, though the demes of the child lineages may differ from that of their parent.Here, a parent of the blue deme sires one child lineage in each of the blue and yellow demes.The production of an event is an integer vector, with one entry for each deme.The production of this event is therefore r " pr blue , r yellow q " p2, 1q.The deme occupancy of an event is the number of lineages in each deme just to the right of the event.The deme occupancy at this event is therefore n " pn blue , n yellow q " p3, 5q.(B) A death-type event causes the extinction of a lineage.Since internal nodes without children are recursively removed, the affected branch is dropped.The production of this event is r " p0, 0q and the deme occupancy is n " p3, 4q.
(C) A migration-type event results in the movement of one or more lineages from one deme to another.Here, one lineage moves from the yellow to the blue deme.The production of this event is r " p1, 0q, i.e., the production is 1 for the blue deme and 0 for the yellow.The deme occupancy is n " p6, 2q.(D) In a sample-type event, one or more sample nodes (blue circles) are inserted.Here, there are two samples, one in each of the blue and yellow demes.Accordingly, r " p1, 1q and n " p2, 6q.(E) A neutral-type event has no effect on the genealogy and zero production in all demes: r " p0, 0q, n " p5, 3q.(F) The theory presented here allows for compound events.As an example, here a birth/death-type event occurs, wherein one yellow lineage is extinguished and a blue lineage simultaneously sires a blue child.For this event, we have r " p2, 0q and n " p6, 2q.(G) Here, a compound sample/death-type event with r " p0, 0q and n " p2, 5q occurs.A blue lineage is sampled and simultaneously extinguished.Note that recursive removal does not occur, since sample nodes are never removed.(H) A compound birth/migration-type event with r " p4, 0q and n " p6, 2q.
(d) At a sample-type event (Fig. 4D), the appropriate number of sampled lineages are selected at random from among the tip nodes, with uniform probability according to deme.One new sample node is introduced for each selected lineage: each is inserted between a selected tip nodes and its ancestor.(e) At a neutral-type event (Fig. 4E), no change is made to the genealogy.(f) Finally, events of compound type (e.g., Fig. 4F-H) are accommodated by combining the foregoing rules.
In each of these events, the new node or nodes that are introduced have node-times equal to the time of the jump.
3.2.1.Emergent lineages and production.The lineages which descend from an inserted node are said to emerge from the event.Thus, after a birth-type event, the emerging lineages include all the new offspring as well as the parent.Likewise, at pure migration-or sample-type events, each migrating or sampled lineage emerges from the event.At pure death-type events, no lineages emerge.In general, at an event of mark u, there are r u i emergent lineages in deme i.We require that r u i be a constant, for each u and i.Thus there is a function r : U ˆD Ñ Z `, such that r u i lineages of deme i emerge from each event of mark u.Since, in applications, one is free to expand the set of jump-marks U as needed, this is not a restriction on the models that the theory can accommodate.We say r u -pr u i q iPD is the production of an event of mark u.Note that the lineages that die as a result of an event do not count in the production but that a parent lineage that survives the event does count.

Conditional independence and exchangeability.
Application of these rules at each jump of X t constructs a chain of genealogies Ĝk .In particular, at each jump-time Tk , the genealogy Ĝk´1 is modified according to the jump-mark Ûk to yield Ĝk .We view Ĝk as the embedded chain of the continuous-time genealogy process G t .It is very important to note that, conditional on p Xk , Ûk q, the number of parents and number of offspring in each deme is determined and the random choice of which lineages die, migrate, are sampled, or sire offspring is independent of these choices at any other times and independent of p Xj , Ûj q for all j ‰ k.Moreover, by assumption, the lineages within each deme are exchangeable: any lineage within a deme is as likely as any other lineage in that deme to be selected as a parent or for death, sampling, or migration.Finally, note that G t does not have the Markov property, though pX t , U t , G t q and pX t , G t q do.3.3.Pruned and obscured genealogies.The process just described yields a genealogy that relates all extant members of the population, and all samples.Moreover, it details each lineage's complete history of movement through the various demes.However, the data we ultimately wish to analyze will be based only on samples.Nor, in general, will the histories of deme occupancy be observable.A generative model must account for this loss of information.We therefore now describe how genealogies are pruned to yield sample-only genealogies and then obscured via the erasure of color from their branches (Fig. 5).

Pruned genealogy.
Given a genealogy G, one obtains the pruned genealogy, P " prunepGq by first dropping every tip node and then recursively dropping every childless internal node (Fig. 5A-B).In a pruned genealogy only internal and sample nodes remain, and sample nodes are found at all of the leaves and possibly some of the interior nodes of the genealogy.Observe that a pruned genealogy is a colored genealogy: it retains information about where among the demes each of its lineages was through time (Fig. 5B).Note also that a pruned genealogy P is characterized by its time, tpP q and the functions P Y and P Z just as an unpruned genealogy is.Finally, observe that, since it contains within itself all of its past history, the pruned genealogy process P t " prunepG t q is Markov, even though the unpruned genealogy process, G t , is not.

Lineage count and saturation.
In the following, we will find that we need to count the deme-specific numbers of lineages present in a given pruned genealogy at a given time.Accordingly, suppose P " pT, Z, Y q is a pruned genealogy and suppose t P r0, T s.Let ℓ i denote the number of lineages in deme i at time t and ℓ -pℓ i q iPD P Z D `.Clearly, ℓ depends only Y t .Therefore, we can define ℓ as a function such that, whenever P " pT, Z, Y q is a pruned genealogy, ℓpY t q is the vector of deme-specific lineage counts at time t.We refer to ℓ as the lineage-count function (cf.Fig. 6).
We will also have occasion to refer to the deme-specific number of lineages emerging from a given event.In particular, given a node time t in a pruned genealogy P " pT, Y, Zq, the number s i of lineages of deme i emerging from all nodes with time t is well defined and we can write sps i q iPD .Like the lineage-count function, s depends only on the local structure of P.However, s depends not only on Y t , but also on r Y t .Thus, we can define the saturation function such that, whenever P " pT, Y, Zq is a pruned genealogy, sp r Y t , Y t q is the integer vector of deme-specific numbers of emerging lineages at time t.Fig. 6 illustrates.

Compatibility.
Suppose P is a pruned genealogy, with tpP q " T and t P evpP q.The local structure of P at t is, in general, compatible with only a subset of the possible jumps U.For example, if the event in P at t is a branch node or a sample node, then it is compatible only with birth-type or sample-type jumps, respectively.Similarly, if the node in P at time t is one at which a lineage moves from deme i to deme i 1 , then u must be either of i Ñ i 1 migration type or of a birth type with parent in i and r u i 1 ą 0. To succinctly accommodate all possibilities, let us introduce the indicator function Q such that Q " 1 if the local genealogy structure-which is captured by the values of P Y just before and after t-is compatible with an event of type u and Q " 0 otherwise.That is, Q u py, y 1 q " 1 if and only if there is a feasible genealogy, G " pT, Z, Yq, and history, H, and a t P r0, Ts such that, given G T " G and H T " H, we have U t " u, r Y t " y, and Y t " y 1 .We refer to Q as the compatibility indicator.The genealogy is pruned by deleting all tip nodes and then recursively pruning away childless internal nodes.Sample nodes are never removed.(C) A pruned genealogy is obscured by effacing all deme information from lineage histories: the colors are erased, as are all inline nodes.See the text ( § §2.9, 3.3.1,and 3.3.4)for more detail.

Obscured genealogy.
The obscured genealogy is obtained by discarding all information about demes and events not visible from the topology of the tree alone (Fig. 5B-C).In particular, if P " pT, Z, Y q is a pruned genealogy, we write obspP q " pT, Zq to denote the obscured genealogy.

Results.
4.1.Likelihood for pruned genealogies.Our first result will be an expression for the likelihood of a given pruned genealogy given the history of the population process.
Theorem 1. Suppose P " pT, Z, Yq is a given pruned genealogy.Define (6) ϕ u px, y, y 1 q -¨npxq ℓpy 1 q r u spy, y 1 q . Lineage count and saturation.Each panel shows the neighborhood of a single event in the unpruned genealogy (top row) and the corresponding pruned genealogy (bottom row).Pruning consists of the removal of all branches that are not ancestral to some sample.
In the bottom row of panels, pruned branches are indicated using broken lines.(A) A birthtype event with production r " pr blue , r yellow q " p1, 1q occurs.(B) Suppose that pruning results in the removal of the dashed lineages.Then the lineage count at this event-time is ℓ " pℓ blue , ℓ yellow q " p2, 2q.The saturation is s " p0, 1q since only a single, yellow lineage emerges from the event.(C) A migration-type event with production r " p0, 1q occurs.(D) After pruning, ℓ " p2, 2q and s " p0, 1q.(E) A sample-type event occurs in which two blue lineages are sampled (production r " p2, 0q).(F) After pruning, ℓ " p2, 2q and s " p1, 0q.
Observe that in panels B and D, the local structures of the pruned genealogies are identical, though they arise from events of different type.
where n is the deme occupancy ( §2.6), r u is the production ( §3.2.1), ℓ and s are the lineage-count and saturation functions, respectively ( §3.3.2),Q is the compatibility indicator ( §3.3.3), and the binomial ratio is as defined in §2.10.Then Prob rP T " P | H T " Hs " 1tevpHq Ě evpPqu Proof.If evpHq ⊉ evpPq, then H and P are incompatible and Prob rP T " P | H T " Hs " 0. Similarly, if any event of H is incompatible with the local structure of P in the sense of §3.3.3, then Prob rP T " P | H T " Hs " 0.
Let us therefore suppose that neither of these conditions hold.Conditional on H T " H, at each time t P evpHq, a jump of mark U t occurred, with a production of r Ut " pr i q iPD , resulting in a deme-occupancy of npX t q " pn i q iPD .In P, at time t, there are ℓ i " ℓ i pY t q lineages in deme i, of which s i " s i p r Y t , Y t q are emergent.By assumption, at each genealogical event, lineages within a deme are exchangeable: each has an identical probability of being involved.This exchangeability implies that each lineage present in a deme at time t was equally likely to have been one of the emergent lineages.In particular, at time t, the probability that s i of the ℓ i deme-i lineages were among the r i of n i lineages emergent in the unpruned genealogy process is the same as the probability that, upon drawing ℓ i balls without replacement from an urn containing r i red balls and n i ´ri black balls, exactly s i of the drawn balls are red, namely ¨ni ´ℓi Because our lineages are labeled, each of the `ℓi si ˘equally probable sets of s i lineages is distinct; just one of these is the one present in P.Moreover, since, again conditional on H T " H, the identities of the lineages involved in a genealogical event are random and independent of the identities selected at all other events, we have established that Prob rP T " P | H T " Hs " Returning to the possibility that H is incompatible with P, since Prob rP T " Ps " 0 if either any Q Ut " 0 or evpPq ⊈ evpHq, we obtain the result.□ Next, we show how the likelihood of a pruned genealogies, unconditional on the history, can be computed.For this, we use the filter equation technology developed in Appendix A. In particular, the following theorem follows immediately from Lemma A2.
Theorem 2. Suppose that P " pT, Z, Yq is a given pruned genealogy.Suppose that w " wpt, xq satisfies the initial condition wp0, xq " p 0 pxq and the filter equation where ϕ is defined in Eq. 6.Then the likelihood of P is LpPq " ż wpT, xq dx.
4.2.Likelihood for obscured genealogies.Our next result concerns the likelihood of a given obscured genealogy conditional on the history.
Theorem 3. Suppose that pT, Zq is a given obscured genealogy.Let q and π be probability kernels, such that for all x P X and y P Y 0 pZq, qpx, yq ě 0, and, for all u P U, t P R `, x, x 1 P X, y, y 1 P Y t pZq, π u pt, x, x 1 , y, y 1 q ě 0, ÿ y 1 PYtpZq π u pt, x, x 1 , y, y 1 q " 1.
Suppose moreover that π u pt, x, x 1 , y, y 1 q ą 0 whenever α u pt, x, x 1 q Q u py, y 1 q ą 0 and that qpx, yq ą 0 whenever Prob " P Y 0 " y ˇˇX 0 " x ‰ ą 0. Then there is a stochastic jump process y t with sample paths in YpZq such that pX t , U t , y t q is Markov and where ϕ is defined in Eq. 6 and the expectation is taken over the sample paths of y t .
Proof.First, observe that, since obs is a deterministic operator, (8) Our strategy will be to evaluate Eq. 8 using importance sampling: we will propose pruned genealogies compatible with Z as sample paths from a stochastic process driven by X t and evaluate the the expectation in Eq. 8 by summing over these paths.Conditional on H T " H, the initial distribution q and probability kernel π generate a Markov chain, ŷk such that Prob rŷ 0 | H T " Hs " qpX 0 , ŷ0 q, Prob rŷ k | ŷk´1 , H T " Hs " π Ûk p Tk , Xk´1 , Xk , ŷk´1 , ŷk q.
The required process y t is the unique càdlàg process with event times Tk and ŷk as its embedded chain.This construction of y t obviously guarantees that evpHq Ě evpyq Ě evpZq and that pX t , U t , y t q is Markov.Now, for y P YpZq, let us define Cpyq " pT, Z, yq.Then, by construction, obspCpyqq " pT, Zq and, conversely, for every pruned genealogy P satisfying tpPq " T and P Z " Z, CpP Y q " P.Moreover, the conditions on the kernels q and π guarantee that, if Prob rP T " P | H T " Hs ą 0 and P Z " Z, then Prob " y " P Y ˇˇH T " H ‰ ą 0. We therefore have that the expectation being taken with respect to the random process y.Here, by definition, πpy|Hq " qpX 0 , y 0 q ź tPevpHq π Ut pt, r X t , X t , r y t , y t q.
The result then follows from Theorem 1. □ Note that, since Y t pZq is finite, it is permissible, for example, to choose q and π to be uniform.
The final result shows how to compute the likelihood of an obscured genealogy.It is an immediate consequence of Theorem 3 and Lemma A2.
Theorem 4. Let V " pT, Zq be a given obscured genealogy.Then there are probability kernels q and π as in Theorem 3 such that if β u pt, x, x 1 , y, y 1 q " α u pt, x, x 1 q π u pt, x, x 1 , y, y 1 q, Ψ u pt, x, x 1 , y, y 1 q " ϕ u px 1 , y, y 1 q π u pt, x, x 1 , y, y 1 q , and if w " wpt, x, yq satisfies the initial condition wp0, x, yq " p 0 pxq 1tqpx, yq ą 0u and the filter equation then the likelihood of V is LpV q " ÿ y ż wpT, x, yq dx.
Lemma A3 shows how this can be computed via Sequential Monte Carlo.

Discussion.
The theory presented here represents a strict generalization of the existing coalescent and birth-death process approaches to phylodynamic inference.In Appendix B, we demonstrate that both of the latter processes are special cases of the genealogical processes constructed here.Importantly, because the theory allows computation of the likelihood via strictly forward-in-time computations, it permits consideration of models for which time-reversal arguments are not available.Moreover, inasmuch as the formulae of Theorem 4 can be efficiently computed via sequential Monte Carlo, explicit expressions for transition probabilities are not needed: it is sufficient to be able to simulate from the population process.This feature of the algorithms-known as the plug-and-play property (He et al., 2010)-further expands the class of population models that can be confronted with data.
In particular, the theory gives us the freedom to choose models with many demes.For deterministic population models, Volz (2012) and Rasmussen et al. (2014) showed how one could accommodate discrete population structure.Their procedures involve solving a large number of differential equations backward in time, relying on the time-reversibility of deterministic dynamics.In general, this time-reversibility is not a property of stochastic processes.
Some existing methods put rather severe limits on the form of the sampling model and, as Volz & Frost (2014) pointed out, misspecification of the sampling model can lead to large inferential biases.With the theory presented here, essentially arbitrary specification of the sampling model is possible.In particular, one can posit sampling at a rate which is an arbitrary function of time and state and include discrete sampling events as well.It is also possible to condition on the existence of samples.
If Sequential Monte Carlo algorithms are used to compute the likelihoods of Theorem 4, then it is straightforward to simultaneously assimilate information from both time-series and genealogical data.One can therefore supplement traditional incidence, disease, or mortality time series with genealogical data in an inferential exercise.
A limitation of the theory is that the population models are assumed to be pure jump processes, which allows consideration of demographic stochasticity and environmental stochasticity modeled by jumps involving multiple individuals (Bretó & Ionides, 2011), but disallows stochastic processes with a diffusive component.It should be possible to incorporate of the full range of Markovian environmental stochasticity via extension of this theory to population models containing both diffusion and jump components.
The price of the theory's flexibility is primarily computational.When Sequential Monte Carlo is used to evaluate the likelihood in Theorem 4, the computational effort scales linearly with the number of samples.In its most straightforward implementation-using an event-driven algorithm (e.g., Gillespie, 1977)-it scales nonlinearly with population size in general.However, stochastic simulation schemes are available that scale independently of population size (Higham, 2008).On the other hand, the importance sampling underlying Theorem 4 will in general require effort that is exponential in the number of demes.For models with many demes, therefore, approaches for ameliorating or circumventing this curse of dimensionality may be necessary.Critically, the substantial freedom one has in the choice of the importance-sampling distribution π can be exploited for this purpose.In particular, since it is permissible to "borrow information" from the future by means of the importance sampling, there is hope for highly efficient algorithmic computation.
In the limit, we obtain Eq.A4, the regular filter equation generated by α, with boost B and zero decay.□ When events are known to have occurred at particular times, it is of interest to integrate over those histories that include an event at each of these times.This leads to singular filter equations, as the next lemma shows.Before we state the lemma, some terminology is needed.Let S be the space of increasing, locally finite sequences in R `, with the topology induced by the Skorokhod metric and Lebesgue measure.For t P R `and s P S, let s t -s X r0, ts.Thus if s P S and s t " pŝ 1 , . . ., ŝK q, then the infinitesimal element of Lebesgue measure at s t is ds t " ś K n"1 dŝ n .Lemma A2.Suppose that B : R `ˆX 2 Ñ R `is measurable and V t is an R `-valued random process satisfying Bpe, r X e , X e q.
Let λ t be a family measures on X ˆS defined by Proof.The proof proceeds by induction on the cardinality of s t .The base case, for which s t " ∅, follows immediately from Lemma A1.Assuming that it holds for |s t | ă K, one has only to verify Eq.A6.This can be accomplished by integrating Eq. 5 directly.□ Remark.In the same way that Eqs. 3 and 4 can be represented as a single equation by means of a Dirac delta notation, the Eqs.A5 and A6 can be collapsed into a more compact form if β is allowed to have atoms at a countable set of time-points and the boost B is adjusted appropriately.
Filter equations afford a convenient means of computing expectations and likelihoods for pure jump processes.This is facilitated by the following Lemma, the statement of which uses a one-sided Dirac delta function.Specifically, let δpv, v 1 q be the right-sided Dirac delta function satisfying δpv, v 1 q " 0 for v ‰ v 1 and ż b a f pvq δpv, v 1 q dv " f pv 1 q 1tv 1 P ra, bqu, whenever f is càdlàg and ´8 ď a ă b ď 8.
Proof.For each t R S, we have Here, the non-explosivity assumption guarantees that we can differentiate under the integral sign and exchange the order of integration.Moreover, it ensures that u Ñ 0 as v Ñ 8. Hence, by evaluating the first integral with respect to v, the second with respect to v 1 , and the third by parts, we obtain which is simplified to obtain Eq.A2.Similarly, at each t P S, we have wpt, xq dx " which is equivalent to Eq. A3 □ Remark.Eqs.A7 are recognizable as the KFE of a certain process pX t , V t q.In particular, the driver X t has KFE Eq.A1.V t is directed by X t in the sense that V has jumps wherever X does: when X jumps at time t from x to x 1 , V jumps by the multiplicative factor Bpt, x, x 1 q ě 0. Between jumps, V t decays deterministically and exponentially at rate λpt, xq.At the known times in S, X jumps according to the probability kernel π and, V jumps by the factor Apt, xq Bpt, x, x 1 q.If we view V t as a weight, then Lemma A3 tells us how the V t -weighted average of X t evolves in time: this average is simply ş wpt, xq dx.Thus, Lemma A3 shows how to integrate Eqs.A5 and A6 in the Monte Carlo sense.Likelihood computation for the SIRS model by Sequential Monte Carlo.(A) A simulated genealogy for β " 4, γ " 2, ω " 1, ψ " 1, pS 0 , I 0 , R 0 q " p97, 3, 0q.(B) A slice through the likelihood surface at the true parameters in the γ-direction.Each point is a distinct Monte Carlo estimate.The blue curve is a LOESS smooth; the dashed lines bound the Monte Carlo-adjusted 95% confidence interval (Ionides et al., 2017).
where N " S `E `I `R is the total population size.Note that the terms associated with sampling cancel each other in the KFE, since, in this model, sampling has no effect on the state.This model has two demes: D " tE, Iu.Its deme occupancy function is npxq " pE, Iq.There are five kinds of jumps: transmission, progression, recovery, waning of immunity, and sampling.The corresponding marks are U " tTrans, Prog, Recov, Wane, Sampleu.Table B2 gives α u , r u , and the event type for each of these marks.

TABLE B2.
Elements of the SEIRS model pertinent to the genealogy process.The table shows the rate (α u ), jump (x Þ Ñ x 1 ), production (r u ), and event type for each of the model's five marks (u).
The filter equation corresponding to the scheme of Table B3 is presented in Box 1.Some numerical results are presented in Fig. B2.TABLE B3.Elements of a scheme for numerically computing the likelihood under the SEIRS model.The regular portion of the filter equation holds in between genealogical events (i.e., for t R B Y S 0 Y S 1 ); the singular portion describes the effect of these events (t P B Y S 0 Y S 1 ).
Each line corresponds to a potential event, but only those for which Q " 1 appear in the equation.An event with mark u and saturation s has the boost given by the binomial ratio shown (third column).The y Þ Ñ y 1 column depicts the proposed painting schematically, and π is the probability of that proposal.Blue is used for the E deme and yellow for the I deme.For each line, the filter equation contains m terms.An asterisk (˚) stands for cases not explicitly mentioned.
Box 1: Filter equation for the SEIRS model As previously, given an obscured genealogy, let B be the set of its branch times, S 0 be the set of tip-sample times, and S 1 be the set of inline sample times.Then for t R B Y S 0 Y S 1 , the filter equation reads: Here, r y k j refers to the coloring of the tree immediately preceding the proposal indicated on line j of Table B3.The integer k specifies the particular branch on which the change occurs.The singular portion of the filter equation has one component for each distinct type of genealogical event: ψ pI ´ℓI q r wpt, S, E, I, R, r y15q, t P S0, ψ r wpt, S, E, I, R, r y18q, t P S1.
In addition to Eqs.B8 and B9, the quantity w should satisfy the condition wpt, S, E, I, R, yq " 0 whenever E ă ℓ E or I ă ℓ I .
A variety of importance-sampling kernels are permissible under the terms of Theorem 4. With a particular choice of importance-sampling kernel, the filter equation uniquely specifies a Sequential Monte Carlo algorithm for estimating the likelihood.The choices made in Table B3 underlie the results displayed in Fig. B2.
B.3.Two-strain competition model.A simple model for the competition of two strains for susceptible hosts is depicted in Fig. 1B.This example will be included in a forthcoming draft.
B.4. Superspreading model.Fig. 1D depicts a model of superspreading.This example will be included in a forthcoming draft.
B.5. Linear birth-death model.In this model, the state variable is the size, N t , of a population at time t.All individuals face the same per-capita birth and death rates, which are λ and µ, respectively.The KFE is Bv Bt " λ pn ´1q vpt, n ´1q ´λ n vpt, nq `µ pn `1q vpt, n `1q ´µ n vpt, nq Stadler (2010) considered the case where samples are taken through time at a uniform per-capita rate ψ.In this case, since there is only one deme, in the filter equation, w can be taken to be independent of y.If B is the set of branch-times and S 0 , S 1 are the sets of terminal and inline samples, respectively, then the regular part of the filter equation is Likelihood computation for the SEIRS model by Sequential Monte Carlo, using the scheme of Box 1. (A) Simulated genealogy for β " 4, σ " 1, γ " 1, ψ " 1, ω " 1, pS 0 , E 0 , I 0 , R 0 q " p200, 3, 5, 100q.(B) Likelihood slice in the σ-direction.Each point represents the estimate of an independent Sequential Monte Carlo computation.The blue curve shows a LOESS smooth; the dashed vertical lines enclose the Monte Carlo-adjusted 95% confidence interval (Ionides et al., 2017).and the singular part is Eqs.B10 and B11 are supplemented by the ancillary condition wpt, nq " 0 for n ă ℓptq.In the Moran model, events occur according to a rate-µ Poisson process.At each event, a compound birth-death jump (cf.Fig. 4F) occurs so that the population size, n, remains constant.If we let X t be the number of events that have occurred by time t, then X t is a simple counting process, which we can use to define the state of the population process.Its KFE is then Bv Bt " µ px ´1q vpt, x ´1q ´µ x vpt, xq, vp0, xq " # 1, x " 0, 0, x ą 0.
Since there is only a single deme, and since nothing depends on the state, in writing the corresponding filter equation, we can take w to be independent of both x and y.
In the classical case (Kingman, 1982a), m samples are taken simultaneously at a single time, T .Then, if B is the set of branch-times and ℓptq is the number of lineages in the genealogy at time t, the filter equation reads   ) shows the genealogy, simulated for λ " 1.2, µ " 0.8, ψ " 1, N 0 " 5. Panel (B) shows a likelihood slice, through the true parameters in the µ direction.As computational effort (i.e., number of particles) increases, the Monte Carlo estimates converge on the exact values, for which an explicit formula exists in this case.The dashed horizontal lines show the approximate maximized likelihood and the 95% critical value (under the likelihood-ratio test).Panels (C-E) are log-log plots that show how the root-mean-square error (RMSE), imprecision (SD), and bias decrease with effort.Note that the bias and the SD are roughly inversely proportional to the effort and its square-root, respectively, as expected.
Integrating Eqs.B12 and taking logarithms yields (B13) log wpT q " k log µ `n 2 ˘´µ `n 2 ˘m ÿ i"m´k where k " |B| is the number of branch-points in r0, T s and the s i -ş 1tℓptq " iu dt are the durations of the coalescent intervals, i.e., intervals between successive branch-points.We recognize Eq.B13 as the expression for the Kingman (1982a) coalescent (e.g., Wakeley, 2009).More generally, if in addition samples are taken according to a rate-ν Poisson process such that the set of sampletimes in the genealogy is S " S 0 Y S 1 , where S 0 , S 1 are the sets of times of terminal and inline samples, Fig. 1D depicts a model of superspreading.There are three demes (D " tE, I L , I H u).

CFIGURE 5 .
FIGURE 5. Unpruned, pruned, and obscured genealogies from a single realization of the genealogy process induced by the SEIRS model depicted in Figs.1 and 2. (A) A realization of the unpruned genealogy process G t is shown at t " 10.Tip nodes, corresponding to lineages alive at time t " 10 are indicated with black points.Blue points represent samples; green points, internal nodes.Branches are colored according to the deme in which the corresponding lineage resided at that point in time: blue denotes E and yellow, I. (B) The genealogy is pruned by deleting all tip nodes and then recursively pruning away childless internal nodes.Sample nodes are never removed.(C) A pruned genealogy is obscured by effacing all deme information from lineage histories: the colors are erased, as are all inline nodes.See the text ( § §2.9, 3.3.1,and 3.3.4)for more detail.

u α u r u
Event type Birth λN 2 pure birth Death γN 0 pure death Sample ψN 1 pure sample B.6.Moran model and the Kingman coalescent.
FIGURE B3.Likelihood computation for the constant-parameter, linear birth-death-sampling model, according to Theorem 4 via Sequential Monte Carlo.Panel (A) shows the genealogy, simulated for λ " 1.2, µ " 0.8, ψ " 1, N 0 " 5. Panel (B) shows a likelihood slice, through the true parameters in the µ direction.As computational effort (i.e., number of particles) increases, the Monte Carlo estimates converge on the exact values, for which an explicit formula exists in this case.The dashed horizontal lines show the approximate maximized likelihood and the 95% critical value (under the likelihood-ratio test).Panels (C-E) are log-log plots that show how the root-mean-square error (RMSE), imprecision (SD), and bias decrease with effort.Note that the bias and the SD are roughly inversely proportional to the effort and its square-root, respectively, as expected.
Pathogen lineages are to be found in hosts within the E and I compartments only.Accordingly, there are two demes: D " tE, Iu.If there is exactly one lineage per host, then the occupancy, npX t q " pn E pX t q, n I pX t qq, is the integer 2-vector giving the numbers of hosts in the respective compartments.
Examples of discretely-structured population models.Demes are shaded.Compartments containing infectious hosts are outlined in green.Curved green lines connect transmission rates with the compartments whose occupancies control their modulation; each such connection gives rise to a nonlinearity in the model.(A) An SEIRS model.Susceptible individuals (S), once infected, enter a transient incubation phase (E) before they become infectious (I).Upon recovery (R), individuals experience immunity from reinfection.If this immunity wanes, they re-enter the susceptible compartment.
λ t pE, Sq " E rV t ¨1tX t P Eu ¨1tDs P S s.t.evpH t q Ě s t us , whenever E Ď X and S Ď S are measurable.Let wpt, x, sq be the density of this measure, i.e., λ t pdx dsq " wpt, x, sq dx ds t ., sq αpt, x 1 , xq Bpt, x 1 , xq dx 1 ´ż wpt, x, sq αpt, x, x 1 q dx 1 ,

TABLE B4 .
Elements of the linear birth-death-sampling model pertinent to the genealogy process.